function [psi] = psi_compute_t(T,q,itrend);
% Compute PSI Matrix for Cosine Transformation
% Inputs"
% T = time series length
% q = number of cosine transforms
% itrend = 1;  Demean
% itredn = 2;  Linear trend
% Output:
% psi = Txq matrix
% Note that here psi'*psi -> I(q)  (dct divided by sqrt(T))

% -- Constants for detrended psi
avec = [
8.986818915818128
15.450503673875414
21.808243318857798
28.132387825662946
34.44151054386154
40.74260591857512
47.03890499737801
53.33210851762535
59.62319758178592
65.91277807964495
72.20124448875121
78.48886472232839
84.77582713626384
91.06226802798255
97.34828846390877
103.63396497455933
109.91935657577787
116.20450950899118
122.4894605207488
128.77423918111484
135.05886955428824
141.343371423239
147.6277612013613
153.91205262066237
160.19625725789024
166.4803849414468
172.76444406945743
179.04844186083437
185.33238455524568
191.6162775737234
197.90012564866376
204.18393292981528
210.4677030712752
216.75143930334937
223.03514449226202
229.31882119004615
235.60247167644877
241.886097994303
248.1697019795255
254.45328528666866
260.7368494107768
267.02039570615676
273.303925402558
279.5874396191708
285.87093937677855
292.1544256083407
298.43789916823965
304.7213608403825
311.00481134532134
317.28825134652686
323.57168145593135
329.85510223883756
336.1385142182776
342.42191787889124
348.70531367038535
354.9887020106256
361.2720832884058
367.5554578659326
373.83882608106046
380.12218824930414
386.40554466565516
392.688895606224
398.97224132972667
405.25558207883364
411.53891808139457
417.8222495515526
424.10557669076024
430.3888996887052
436.6722187241573
442.95553396574263
449.23884557265296
455.5221536952966
461.8054584758957
468.088760049036
474.3720585421723
480.6553540760948
486.93864676535895
493.2219367186826
499.5052240393135
505.78850882536835
512.0717911701483
518.3550711624309
524.63834888674
530.9216244235986
537.2048978497622
543.4881692384369
549.7714386594818
556.0547061795986
562.3379718625079
568.6212357691134
574.9044979576566
581.1877584838604
587.4710174010634
593.7542747603464
600.0375306106511
606.32078499889
612.6040379700509
618.8872895672948
625.1705398320477
631.4537888040863
];

  fvec = pi*(1:1:q);
  tvec = (0.5:1:T)'/T;
  dct = (sqrt(2)/sqrt(T))*cos(tvec*fvec);
  fvec_2T = fvec/(2*T);
  ljt = (sin(fvec_2T))./fvec_2T;
  
  if itrend == 1;
    psi = dct.*repmat(ljt,T,1);
  elseif itrend == 2;
    % Compute Demeaned Version of DCT
    z = [ones(T,1) (1:1:T)'];
    zzi=inv(z'*z);
    b=zzi*(z'*dct);
    dct_t=dct-z*b;
    % Compute asymptotic eigenfunctions for detrended random walk -- discarding first two 
    ef_t=dct_t;               % Note ... even number columns of ef_t are the same as dct_t, same as dct 
    mjt=ljt;
    for i = 3:2:q;
      sgn=(-1).^((i+1)/2);
      anum=(i-1)/2;
      om=avec(anum);
      ef_t(:,i)=sgn*sqrt(2*om/(om-sin(om)))*(sin((tvec-0.5*ones(T,1))*om)/sqrt(T));
      tmp=om/(2*T);
      mjt(i)=sin(tmp)/tmp;
    end;
    b=zzi*(z'*ef_t);
    ef_t=ef_t-z*b;               % This imposes exact orthogonality for this value of T 
    ef_t=ef_t(:,2:size(ef_t,2)); % Drop first column of ef_t and dct_t ... start at frequency 2*pi/t @
    mjt=mjt(1,2:size(mjt,2));
    psi=ef_t.*repmat(mjt,T,1);
  end;
    
end
